Kinetics of a single trapped ion in an ultracold buffer gas 
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The immersion of a single ion confined by a radiofrequency trap in an ultracold atomic gas extends 
the concept of buffer gas cooling to a new temperature regime. The steady state energy distribution 
of the ion is determined by its kinetics in the radiofrequency field rather than the temperature of 
the buffer gas. Moreover, the finite size of the ultracold gas facilitates the observation of back- 
action of the ion onto the buffer gas. We numerically investigate the system's properties depending 
on atom-ion mass ratio, trap geometry, differential cross-section, and non-uniform neutral atom 
density distribution. Experimental results are well reproduced by our model considering only elastic 
collisions. We identify excess micromotion to set the typical scale for the ion energy statistics and 
explore the applicability of the mobility collision cross-section to the ultracold regime. 

PACS numbers: 34.10.-l-x 37.10.Ty 34.50.-s 05.10.Ln 



I. INTRODUCTION 

Trapped ion systems are among the most promising 
candidates for quantum information processing [1, 2], 
precision measurements [3], and quantum chemistry [4]. 
For many of these applications it is required to cool the 
ions to low temperatures. To this end, various techniques 
such as laser cooling, resistive cooling, sympathetic cool- 
ing by other ions, or buffer gas cooling are routinely used. 
Ultracold atomic gases have recently become available in 
hybrid systems with trapped ions [5-10], extending the 
concept of buffer gas cooling to ultralow temperatures. 
Understanding this new regime and how it relates to con- 
ventional buffer gas cooling is essential for any future 
application to trapped ions. 

Cooling of an ion in a buffer gas is caused by elastic 
collisions. They are dominated by the long range polar- 
ization interaction, and cross-sections are considerably 
larger than for collisions between two neutral atoms [11- 
21]. Therefore, collision rates are large and cooling is 
expected to be efficient [22]. Moreover, it has been pro- 
posed that internal degrees of freedom of ionic molecules 
could also be cooled by using an ultracold buffer gas [23] . 
Another specific feature of the polarization interaction is 
that collisions affecting the ion's mobility happen with 
rates independent of the collision energy [24]. This can 
lead to simplified system behaviour and has, for example, 
been applied in ion mobility spectrometry [25]. 

The motion of an ion in a radiofrequency (RF) trap 
can be decomposed into a fast driven motion, the mi- 
cromotion, and a slow secular motion. In every collision 
the energy of the RF-field couples via the micromotion 
to the neutral atom's energy and the ion's secular energy. 
This can lead to ion energy removal or intake, sensitive 
to the RF-phase in the moment of the collision [26-28]. 
The average effect of many collisions results in cooling 
or heating and depends on parameters such as the ratio 
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between the mass of the ion and the mass of the neu- 
tral atom m„. For very heavy neutrals runaway heat- 
ing of the ion is expected [26], whereas very light neu- 
trals enable efficient buffer gas cooling. In the mass ra- 
tio regime between the two extremes increased ion trap 
loss can be observed [29, 30]. This effect has been ex- 
plained by non-thermal energy distributions of the ion 
obtained from Monte Carlo simulations [31]. Such nu- 
merical simulations are a well established tool to model 
a trapped ion interacting with a buffer gas [32-34] . They 
can account for energy dependent scattering rates, com- 
plex electric field geometries, and other experimental pa- 
rameters, which are difficult to treat analytically. In pre- 
vious calculations, buffer gases at ambient temperatures 
have been assumed, and the cooling of the ion has been 
limited to the buffer gas temperature. However, in the 
recent experiments with ultracold neutral buffer gases a 
new energy scale related to excess micromotion has be- 
come dominant. The direct relation between the ion's 
mean energy and the excess micromotion has been ob- 
served in [8], for a system of Yb"*" — Rb. 



Here, we investigate the kinetics of a single ion collid- 
ing elastically with an ultracold buffer gas, by applying 
Monte Carlo techniques. The effects on the ultracold 
neutral cloud are modelled using a semiclassical differen- 
tial cross-section. The results on neutral atom loss and 
temperature increase, and the dependence of ion energy 
on excess micromotion are in good agreement with ex- 
periments. 



The paper is organized as follows: In section II we 
describe the basic simulation procedure and the under- 
lying physical model. The classical Langevin interaction 
model is applied in section HI to derive the ion's energy 
statistics depending on the mass ratio, trap geometry, 
and scattering rate. In section IV we explain effects on 
the neutral atom cloud as a result of the energy depen- 
dent differential cross-section and the non-uniform neu- 
tral density distribution. 
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II. SIMULATION MODEL 

The time evolution of a trapped single ion colliding 

with ultracold atoms is modelled using a simulation con- 
sisting of an analytical and a numerical part. In the time 
between collisions, trajectories are analytically described 
using the pseudo-potential approximation, while clastic 
collisions are taken into account using Monte Carlo tech- 
niques. 



A. Ion trajectory 

We consider a single ion confined by the RF- 
quadrupole potential of a linear Paul trap 



r 



(1) 



Here, Vo is the RF-voltage amplitude applied with fre- 
quency Q,T to electrodes at a distance Rt away from the 
trap symmetry axis. In addition, we consider a static 
quadrupole potential confining the ion along the trap 
symmetry axis with trapping frequency oJz, 



^static — 2 Q ' 
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meaning, the timescale of the collision is shorter than 
ri^^, which is the shortest timescale of the motion of the 
ion. This assumption implies that every collision is sen- 
sitive to the momentary relative velocity, including the 
micromotion [37]. Therefore, we consider a number of ef- 
fects causing contributions to the micromotion. Firstly, 
the intrinsic micromotion described in Eqn. (3) is propor- 
tional to the distance of the ion from the centre of the 
RF-quadrupolc field. Secondly, static offset electric fields 
displace the minimum of the ion trapping potential by a 
distance (Ax, Ay) from the symmetry axis of Eqn. (1). 
Thirdly, RF pickup on end-cap electrodes can lead to 
micromotion along the trap symmetry axis. These con- 
tributions are summed in the expression 



-rsec,y{t) - /^y \cos{Q.Tt), (5) 



where Cz parameterizes the micromotion along the trap 
symmetry axis [38]. 

The ion velocity considered for collisions is 



(6) 



with Vsec = ^'f^sec- This is similar to the time derivative 
of Eqn. (3) but includes all the excess micromotion terms. 



Q is the ion's charge. Mathieu equations describe the 
classical motion of an ion in the combined potential (see 
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for example [35, 36]), using the parameters a = 2 ^ and 



9 = ^/8 ^ with a;, = ^ For a < /2 « 1 

the Floquet solution to first order in q yields 



= Ay Sm{u)yt + iPy) 
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(3b) 



which is usiially referred to as the pseudo-potential ap- 
proximation. It consists of a rapidly oscillating mi- 
cromotion term and the secular motion, which is har- 
monic with frequencies ujx,y = — ^wf and ampli- 

1 I '2, E — ♦ 

tudes Ax,y = ^j^y ^.'^ ■ A full secular trajectory fsec 
including the harmonic motion along the trap symmetry 
axis is described by three energies Ej and three phases 
'Pj, j e {x,y,z}, with 



sec, J 



2E^ 



smiijjj t + ifj) . 



(4) 



This formula, describing a three-dimensional harmonic 

oscillator, will be used throughout the following calcula- 
tions to approximate the ion's position. The total secular 
energy Ex + Ey + E^ will be referred to as the ion energy. 

The motion of the ion is affected by collisions with the 
neutral atoms. We assume them to be instantaneous, 



B. Collision dynamics 

The simulation uses classical trajectories for the mo- 
tion of the ion. Therefore, its validity is restricted to 
ion energies well above the energy quanta of secular mo- 
tion fko. The temperature of the ultracold buffer gas is 
assumed to be well below this energy scale, and in the 
collision the neutral atom's initial energy is neglected. 
Due to conservation of energy and momentum, the elas- 
tic scattering process is defined by the scattering angles 
{0, (j)). The ion's velocity changes according to 



VionJ = (1 - /3) VionS +1311 Vu 



(7) 



with Vion,i (vionj) being the initial (final) velocity given 
by Eqn. (6) at the time of the collision, (3 = and 



TZ is the rotation matrix determined by 9 and (f), with 
respect to the direction of Vion,i- From Vsecj and fsec a 
new set of (pj and Ej can be determined, which describes 
the ion's trajectory after the collision. 

To illustrate its impact on the motion of the ion, 
Eqn. (7) can be rewritten in terms of the secular velocity, 
yielding 



= (1 - /3) Vsec,i +011 Vsec,i + P {U - 1) Vr, 



(8) 



This expression shows how Vmm couples to the secular 
velocity in every collision. Note that Vmm is, by definition 
of Eqn. (5), the same before and after the collision since 
it only depends on the position fsec and time t of the 
instantaneous collision. 
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C. Scattering rate 

The probability dPc for the ion to colhde with a neutral 
atom within a short time interval dt defines the scattering 
rate 

dP 

r(i) = ^ = n{x) a{E,) v,on{t) . (9) 

It is proportional to the neutral atom density n{x) at 
the ion's position. The cross-section a{Ec) is usually a 
function of the collision energy Ec, which, neglecting the 
energy of the neutral atom, is given hy Ec = P ^ ''^lon- 

The ion's position changes on a timcscale of ujj^ while 

the velocity of the ion Vion changes on a timescale of 17^^. 
In general, V{t) will therefore be time dependent in a 
non-trivial way. The sampling method used to efficiently 
choose the time of collision is explained in Appendix A. 

A technical description of the main simulation loop is 
given in Appendix B. 



D. Inelastic collisions 

Inelastic processes like charge exchange, spin exchange 
or molecule formation have been predicted to occur in 
the hybrid system [11, 22, 39-41]. In experiments [6- 
9], charge exchange, which is typically signaled by the 
loss of the ion, has been observed at rates many orders 
of magnitude lower than the elastic collision rate, in the 
non-rcsonant case. Spin exchange collisions can occur 
with rates comparable to elastic scattering [22], and en- 
ergy from internal states can be transferred to the exter- 
nal degrees of freedom. In a spin-stretched configuration, 
however, spin exchange is suppressed. 

In our simulation, we can include inelastic eff'ects by 
introducing additional, competing rates, defined as in 
Eqn. (9), but with inelastic cross-sections ai{Ec) instead. 
The effect of any inelastic event on the hybrid system de- 
pends primarily on the question whether the original ion 
still exists after the process. If this is not the case, the 
simulation can be stopped at the first occurrence. If the 
ion continues to exist, the amount of energy released or 
absorbed by the internal states of the colliding particles 
needs to be considered in a modified version of Eqn. (7). 
In either case, our simulation is able to predict the rate 
at which inelastic events occur, given the inelastic cross- 
section ai{Ec). On the other hand, the simulation can be 
used, if a rate is measured experimentally, to determine 
<Ji{Ec). 

In the following sections we consider clastic processes 
only, assuming inelastic processes either to happen very 
rarely, in line with the experimens [7-9], or to involve 
only small amounts of internal energy, which do not sig- 
nificantly affect the system. 



III. LANGEVIN SCATTERING 

For the motion of an ion in a neutral gas, mainly 

large angle scattering is considered relevant, as small de- 
flections do not signiflcantly change the ion's trajectory. 
This assumption leads to the Langevin scattering model, 
which successfully describes the ion's mobility in previ- 
ous experiments with ions in a neutral buffer gas. Here, 
we apply the Langevin scattering model to the trapped 
ion system including the effects of micromotion. We in- 
vestigate the properties of the ion's energy and will later 
compare these results with those from a more complete 
semiclassical scattering model, in section IV. We will 
indeed find good agreement between the two models in 
describing the energy distribution of the ion, confirming 
the above assumption, that large angle scattering events 
determine the evolution of the ion's energy. 

The ion-neutral interaction is dominated by the attrac- 
tive polarization interaction, which is of the form 

nfl) = -^ (10) 

with d = aoQ"^ / (iireo)^ being proportional to the neu- 
tral particle polarizability ao- i? is the internuclcar sep- 
aration. Classically one can define a critical impact pa- 
rameter be = {2Ci/ EcY/"^ [42]. Collisions with impact 
parameter b > be lead to small defiections and are ne- 
glected. Impact parameters b < be result in inward- 
spiralling trajectories, which lead to almost uniformly 
distributed scattering angles into all directions as in hard- 
sphere scattering. The resulting cross-section for large 
angle scattering, cr^ = nb"^ is proportional to the inverse 
collision velocity, and leads to a scattering rate indepen- 
dent of the collision energy [24]. 

A. Energy scale 

Our aim is to determine the energy scale of the ion 
on its classical trajectory after many collisions such that 
the initial conditions for Ej can be neglected. Neither 
the Langevin scattering at its energy-independent rate, 
nor the ultracold neutral buffer gas, which we assume at 
T = 0, introduce such a scale. As a consequence the 
only energy scale in the system is defined by the excess 
micromotion, see Eqn. (5). This is in contrast to the case 
of a buffer gas with non-negligible temperature, where 
the energy scale is rather set by the temperature of the 
neutral gas [31]. 

To associate the excess-micromotion with an energy 
scale we define 

_ IJh 2 /I 1 \ 

with Vmm,e being the full velocity amplitude of the mi- 
cromotion for an ion in the centre of the ion trapping 

potential. g depends on the displacement (A.t, Ay) 
caused by the uncompensated offset electric field. Note 
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that when the offset electric field is compensated using 
a photon correlation measurement [43, 44], photon shot 
noise usually limits the lowest achievable Emm.e- 

Since Emm,e is the only energy scale in the colliding 
system it is sufficient to express all energies in units of 
Emm.e- This givcs general results for any amplitude of 
uncompensated micromotion. It also implies that any 
statistical measure of ion energy has to scale with Emm,e 
and therefore with the square of the excess micromotion 
amplitude. Such quadratic dependence of the mean ion 
energy has been experimentally observed in [8]. 

B. Energy spectrum 

In order to treat scattering in the presence of micromo- 
tion in its most general form we choose all trap frequen- 
cies {u)x,y,ZT^T) to have irrational ratios. We assume a 
very low homogeneous neutral atom density such that the 
scattering rate T is much smaller than the trap frequen- 
cies. This ensures that two consecutive collisions happen 
at uncorrelated positions. 

We obtain the energy spectrum from the simulation 
by binning the secular energy after each collision on a 
logarithmic scale. This measures a logarithmic energy 
probability distribution dP(i?)/dlog(£'). Fig. 1 shows 
such energy spectra for three different mass ratios of the 
atom and the ion. For heavier neutral atoms the mean 
energy of the ion increases and a tail in the spectrum to- 
wards higher energies becomes dominant. Even for equal 



masses {j3 = 0.5) the tail towards high energies is evident 
when compared to the thermal distribution with the same 
mean energy. For any mass ratio, the obtained spec- 
trum distinctly differs from a thermal distribution, also 
because very low energies {E <^ Emm,e) are extremely 
rare. This supports the validity of the classical treatment 
of trajectories and instantaneous collisions {E > H^It)- 

A power law, dP{E)/dlog{E) oc E" nicely fits the tail 
in the spectrum towards high energies for m„ > or 
/3 > 0.5 [31]. As /3 is increased further, towards even 
heavier neutrals, the negative exponent a approaches 0, 
at which point the spectrum does not converge with time 
anymore and runaway heating starts to dominate the evo- 
lution of the ion's energy. The critical mass ratio param- 
eter Pcrit can be found by extrapolating the exponent 
a{P) towards a{l3crit) — 0. The quantity jicrit is not a 
universal number but is a function of the trap geometry. 
It depends on the ratio or, expressed in the ion trap 
parameters a and g, on 2- — (Hii)2 rpj^g dependency 
is explained by the fact that axial confinement leads to 
radial deconfinement and thereby to an increase of the ra- 
tio between average micromotion and secular energy. For 
three different cases the extrapolation to Pcrit is shown in 
Fig. 2. Our data are compatible with the critical mass ra- 
tio previously found for a specific trap geometry [31]. For 
very elongated traps (cj^, ^ ujp) we find l^crit — 0.685, 
corresponding to a mass ratio — = 2.17 . 




10 " 10° 10' 10° 10' p 



FIG. 1. Energy spectra of the ion for three different mass ra- 
tios. 10* energies are sampled into logarithmically spaced bins 
with an energy resolution of 100 bins per decade. We choose 
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Ax = Ay and Cz = and trap parameters ^ — 50. In the 
case where the ion mass is twice the neutral mass (black) the 
ion's mean energy is 0.8 -Emm, e, in the case for equal masses 
(red) 5Emm,e- The red dashed line corresponds to a thermal 
energy distribution with 5Emm,e- The spectrum for a lighter 
ion, — 1.7 (blue), contains a significant contribution of 
very high energies, typically leading to quick ion loss due to 
finite trap depth. 



FIG. 2. The exponents a obtained from fitting power laws 
to the high energy tails of the ion energy spectra are plotted 
against the mass ratio coefficient /3. This is done for three 

2 

different trap geometries, for a spherical trap (— — 6, black 

2 

rectangles), for an elongated trap (— = 50, red circles) and 

2 

for the extreme case with negUgible axial confinement (^ = 
10^°, blue triangles). The data is fitted with second order 
polynomials to extrapolate to Peru for which the exponent a 
becomes 0. 
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C. Average energy and lifetime of tlie ion 

To evaluate the effectiveness of buffer gas cooling for 
different mass ratios we calculate the average of the en- 
ergy spectrum of the ion. As a physical measure, the 
arithmetic mean comes close to the definition of a tem- 
perature, albeit the clearly non-thermal distribution. We 

2 

show the arithmetic mean in Fig. 3 for the case ^ = 50. 
Although the arithmetic mean diverges for a > — 1, the 
energy spectrum can still be normalized for a < 0. In 
this range the median continues to be a well defined sta- 
tistical measure for the energy of the ion. 

In experiments, a large probability density at high en- 
ergies leads to a rapid ion loss due to a finite trap depth 
Etd [30, 31]. We have numerically evaluated the required 
trap depth Etd to limit the ion loss probability per col- 
lision to Pioss (Fig. 3). For large /3 we can approximate 

the required trap depth by Etd = Emm^e (Pioss)^^"- The 
results show that mass ratios with light neutrals are pre- 
ferred for efficient buffer gas cooling. However, even for 
the heavy neutral scenario stable trapping and buffer gas 
cooling are possible for carefully chosen trap geometry, 
trap depth and micromotion compensation. 
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FIG. 3. For different /3 the ion's arithmetic mean energy 
(black rectangles) and median energy (red circles) are shown 
in units of Smm.e- The arithmetic mean is expected to di- 
verge for P > 0.554. An example for the required trap depth 
for Pioss < 10~^ is given (blue triangles) and compared to 
Etd ~ Emm,e {PioBs)^^°' (dotted line) using a obtained from 
power law fits to the tails of the energy spectra. All data are 
for 2^ = 50. 



D. Higher collision rates 



at correlated positions increases. Under this condition we 
observe a reduced median energy and an increased power 
law tail. 



IV. SEMICLASSICAL SCATTERING 

Up to this point the classical Langevin model has been 
used to explore the ion's energy spectrum and its depen- 
dence on mass ratio, trap geometry and collision rates. 
Here we make use of a semiclassical description of the 
interaction process. The solution to the quantum me- 
chanical scattering problem can be found by expanding 
the wavefunction into partial waves. In the regime where 
many partial waves contribute, the total elastic cross sec- 
tion scales like Ec [39]. The resulting energy depen- 
dent scattering rate and angular dependence lead to addi- 
tional effects as compared to the Langevin model. These 
are necessary to explain the back action on the neutral 
cloud. 

We model the interaction potential by the long range 
polarization interaction of Eqn. (10) plus repulsion at 
short distances. The full differential cross-section is cal- 
culated using [45] 



da 



1 



OO 

E 

1=0 



(2Z-f 1)6^"' sm{r]i)Piicoa0) 



(12) 



The angular momentum of a partial wave is M and 
hk = ^y2fiEc is the collision momentum. The scat- 
tering phase rji can be obtained by solving the radial 
Schrodinger equation which involves the centrifugal po- 



tential 



2^ji2 — • The resulting centrifugal barrier in- 
creases in height with angular momentum (~ {Hl)'^). Par- 
tial waves with I < Iq = 1/ h\/ l^i^flC^Ec have a collision 
energy larger than the height of the centrifugal barrier, 
probe the deep potential well and are refiected from the 
hard core. The exact form of the potential, relevant to 
determine 7y;, is typically not known. Therefore phase 
shifts r\i for / < Iq are assumed to be uniformly dis- 
tributed within [0,27r) [11]. In this approximation each 
partial wave contributes with ai = ^ to the total cross- 
section. Summing ai up to reproduces the Langevin 
cross-section. 

The full quantum mechanical cross-section includes ad- 
ditional contributions from partial waves with I > Iq. As 
the centrifugal barrier is higher than the collision energy, 
these partial waves are scattered from the centrifugal bar- 
rier, if tunnelling effects are neglected. The phase shifts 
can be semiclassically approximated by [11, 39] 



For all the results discussed so far the collision rate 
has been assumed very low, and as long as the condi- 
tion F ^ is fulfilled the results remain unchanged. As 
the collision rate approaches or even exceeds the trap fre- 
quency i-Ox,y, the probability to have consecutive collisions 



with Rq 



1 + 1/2 
k ■ 



V{R) 
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fc2 - 



('+1/2)" 

— w~ 



(13) 
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A. Modeling the differential cross-section 

The probability distribution for the deflection angle 

is numerically calculated using Eqn. (12), Eqn. (13) and 
randomly distributed rji for I < Iq. We sum Eqn. (12) 
for I up to 20000 and average over 100 different random 
sets of rji . The differential cross-section calculated in this 
way depends on the reduced mass fi, the collision energy 
Ec = -^Y^ ^'^d C4 only. For the case of a ^ '''^Yb''" ion col- 
liding with a *^Rb atom, four probability distributions of 
the form Eqn. (14) are shown in Fig. 4. The main feature 
of I{9,Ec) is a forward scattering peak, which gets more 
pronounced as the energy increases [41]. The integral 
of the differential cross-section reproduces the expected 

— 1/3 

Ec energy dependence, and its magnitude is in agree- 
ment with [11]. 




FIG. 4. (a) The probability distributions 1(6, Ec) for the 
scattering angle 9 in elastic collisions between ^'^^Yb^ and 
*^Rb are numerically calculated. The unknown scattering 
phases for close encounters (/ < Iq) are chosen randomly. 
The four curves are for different collision energies Ec- The 
forward scattering peak at small 9 is more pronounced for 
higher energies. Its shape is emphasized in (b) where the 
scattering angle 9 is displayed logarithmically. The normal- 
ized I{9, Ec) (smooth lines) for the four different energies are 
compared with logarithmically binned output from the ran- 
dom ^-generating function (step like lines). The approxima- 
tion is optimized to reproduce the height and position of the 
forward scattering peak. 

To implement the differential cross-section in the 
Monte Carlo simulation a parameterization of the nor- 
malized I{9, Ec) is used to create a function that returns 
a random 9 for a given collison energy. The distribution 
I{9, Ec) is modelled in four intervals using two power laws 
(oc 9P), a flat top and a flat background. The parameters 
peak height, background offset and interval limits are en- 
ergy dependent and are well approximated by power laws 



(oc EP ). These power laws are obtained from fits to dif- 
ferential cross-sections for more than 30 different energies 
Ec in the range between fcs x 1 /jK and ks x 100 K. The 9- 
generating function uses these parameter functions and 
inverse transform sampling. Fig. 4b compares sampled 
output of the ^-generating function with the normalized 
differential cross-sections. 



B. Effects of the energy dependent scattering rate 
on the ion energy spectrum 

We have simulated the kinetics of the ion in an ultra- 
cold buffer gas using the parametrized differential cross- 
section to investigate how this affects the ion energy spec- 
trum. Different from the Langevin case, the collision rate 
depends on the instantaneous ion energy. Therefore ef- 
ficient simulation relies on the collision time sampling 
described in Appendix A. The ion energies are binned 
and weighted by the time the ion remains at the specific 
energy. 

We have performed simulations for an ion trap with 
frequencies uJx.y.z — 2tt x {151, 153, 42} kHz, Ht — 2tt x 
42.5 MHz, excess micromotion parameters Ax — Ay = 
2 /ini and neutral density n = 10^* m~'^ for the sys- 
tem ^^''Yb^ - ^^Rb, reproducing conditions comparable 
to [8]. We find good agreement with the previous re- 
sults from the Langevin model and conclude that the 
Langevin model is sufficient to describe the ion's energy 
statistics. Formally, however, the system does not nec- 
essarily scale only with the excess-micromotion energy 
Emm,e anymore, as the differential cross-section intro- 
duces its own energy scale. In the next section we will 
demonstrate that the full differential cross-section is nec- 
essary to predict effects on the cold neutral atoms. 



C. Neutral cloud evolution 

Ultracold neutral atomic clouds have atom numbers 
ranging up to 10^. Compared to a room temperature 
buffer gas, this limited number of atoms and the good 
isolation from the environment allow the observation of 
collision effects on the neutral gas. The main observables 
are the number of neutral atoms Na and their temper- 
ature Ta- In experiments, these values can be obtained 
from time-of-flight imaging and are suitable to verify the 
simulation model. 

The back-action of the ion onto the neutral gas is a 
result of the energy transfer per collision. For general 
two-body elastic collisions it is given by 

Et = A{l- p)EcSiiv^{9/2) , (15) 

depending on the scattering angle 0. For very small de- 
flections <C 1, resulting from the forward scattering 
peak in the differential cross-section, only very little en- 
ergy is transferred to the neutral atom. The distribution 
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of transferred energies Et, shown in Fig. 5, can be under- 
stood as a convolution of the cohision energy distribution 
and the energy dependent differential cross-section. 




pending on the atoms energy. The new temperature is 
calculated 



Taj 



Ng^i 3 ks Tg.^ - I fc_B Tg^i - Epot jr) 



(16) 



using the total energy of the neutral cloud, and the poten- 
tial and average kinetic energy of the lost atom known 
from the position r of the collision. This can lead to 
evaporative cooling or heating effects depending on the 
position of the ion in the neutral gas. 

For ion trajectories larger than the size of the buffer gas 
the ion can only collide in the centre of the trap where its 
motional energy is mostly related to the secular motion 
rather than the micromotion. This suppresses the power 
law tail of the ion energy distribution and reduces the 
ion's average energy. Hence, using tight traps to confine 
the neutral atoms might help to overcome the constraints 
on ion trap depth, trap geometry and mass ratio. 



FIG. 5. Distributions of collision energy Ec (dotted line) and 
transferred energy Et (solid line). The data are obtained by 
binning 10* collisions on a logarithmic energy scale, compar- 
ing Langevin scattering (red) with the full semiclassical dif- 
ferential cross-section (blue). The vertical axis indicates the 
scattering rate in kHz per decade of energy at which collisions 
with the specific energies occur. The inset is a magnification 
of the region relevant for Langevin scattering. The settings 
for this simulation run are n = 10^* m~^ for the uniform 
neutral atom density and Ax = Ay = 2 ^m, resulting in an 
excess-micromotion energy of _Emm,e/fcs = 160 mK. Note the 
similarity between the two models for large and the significant 
difference for small transferred energies Et- 



D. Comparison to experimental data 




Fig. 5 also shows distributions for E^ and Et obtained 
using the Langevin model for comparison. The Langevin 
collision rate is obviously smaller, explained by the differ- 
ent energy dependence of the cross-sections, ctl oc Ec 

— 1/3 

vs (T (X i?c • The semiclassical distribution of E^ is 
also slightly shifted towards higher energies with respect 
to the Langevin E^-, since collisions are more likely to 
happen at higher energies. The transferred energies dif- 
fer only little between to two models for Et fcs x 0.03 K. 
This reflects the statement that the ion's mobility is well 
described by Langevin type collisions. However, the sig- 
nificant peak at low Et of the semiclassical distribution 
causes most of the effects on the cold neutral gas. 

So far all the results were obtained assuming a uniform 
density distribution of neutral atoms. This will now be 
replaced by a spatial distribution for a thermal gas in a 
harmonic trap with temperature Tg and atom number 
Ng. Considering a finite trap depth Etd.g for the neu- 
tral atoms, every collision with Et > Etd.g will lead to 
an atom loss, whereas every Et < Etd,a will increase the 
temperature Tg. The simplified model used in the sim- 
ulation assumes immediate thermal equilibration of the 
neutral gas. Then a collision with Et < Etd,a simply in- 
creases Tg by Y^^TT- '^^^ ^'^^^ '^^ atom for Et > Etd,g 
will decrease Ng%y 1 but also affect the temperature de- 



0.0 1^ — I — ' — ^ — ^ — ' — ^ — ^ — ' — ^ — ' — ^ — I I 

01234567 01234567 

^mm.e'^B [^l E^m.e ' [K] 

FIG. 6. Comparison between experimental measurements 
(black circles) and the simulation predictions for neutral atom 
loss (a) and the neutral temperature increase (b). The semi- 
classical model (blue) fits the data well. The Langevin model 
(red) systematically underestimates the collision effects on 
the neutral atoms. The contribution of evaporative heating 
(Eqn. (16)) in the semiclassical model is indicated with the 
blue dashed line. The experimental data are taken from [8]. 

Here we compare the simulation predictions for the 
effects on neutral atoms with experimental data from 
Yb"*" — Rb [8] . The measured quantities are the loss of 
neutral atoms and temperature increase of a cold ther- 
mal cloud, after 8 s of interaction and for different excess- 
micromotion energies E„im.e- Fig. 6 displays the data to- 
gether with the simulation results. 

Initial conditions for the neutral ® ''Rb cloud in the sim- 
ulation are given by Ta,o = 250 nK and Ng^o — 2.25 x 10^. 
Neutral trap frequencies of 27r x {28, 28, 8} Hz result in 
an initial central density of n(0) = 1.9 x 10^® m~^. The 
neutral trap depth is Etd.g — kB 'x S fiK. 

A single -'^''^Yb''" ion is trapped with parameters given 
in section IV B. The excess micromotion parameters Ax 
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and Ay are varied between and 15 /im, and addition- 
ally there is micromotion along the trap symmetry axis 
with C2 = 2 fim. The exeess-mieromotion energy scale, 
Emm,e/kB thus rangcs from 90 mK to 7K. 

The semiclassical simulation predicts both the shape 
and the magnitude of the experimental results well. In 
contrast, the Langevin scattering is not suited to model 
these measurements because the ultracold atoms are 
highly sensitive to small energy transfer Et. They cor- 
respond to collisions with small deflection angles, which 
are neglected by the Langevin model, cf. Fig. 5. 



V. CONCLUSION 

We have numerically investigated the kinetics of a sin- 
gle trapped ion interacting with an ultracold neutral gas. 
Our results explain the effects of the mass ratio, trap ge- 
ometry and excess micromotion on the ion's energy spec- 
trum. We have applied two different collision models 
of the atom-ion interaction, the Langevin and the semi- 
classical scattering model. Both yield similar ion energy 
statistics. The greater simplicity of the Langevin model 
introduces a characteristic energy scale for the ion en- 
ergy statistics. However, experimentally observed effects 
on the neutral cloud can only be explained by the semi- 
classical model. Forward scattering events with small 
energy transfer are affecting both the neutral cloud tem- 
perature and the atom loss rate. Cold atom-ion collisions 
could then be used to locally remove atoms resulting in 
efficient cooling of quantum gases. 
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Appendix A: Collision time sampling 

A scattering rate F sets the probability for a collision 
to take place within a time interval dt. In our case 

T{t) = n{x) a{Ec) it) (Al) 

with n being the neutral atom density, a the cross-section 
and the velocity Vion as defined in Eqn. (6). Usually the 
density is a function on the; ion's position, the cross- 
section depends on the collision energy and Vion oscillates 
rapidly in time, leading to a r{t) with non-trivial time 
dependence. In the following we explain the method used 
to randomly generate collision times t with the exact dis- 
tribution defined by T{t). 

In general, the process of an event (collision) taking 
place with a rate T{t) can be modelled using a differential 



equation for the probability Q{t) that the event has not 
yet happened after the time t, 

dQ{t) = -r{t) Q{t) dt . (A2) 

The probability distribution for an event to take place 
after the time t is defined by P{t) = —dQ/dt. In the 
simple case with constant r(t) = Fq the solution is 

Prjt) -Fo exp(-Foi) (A3) 

and a random time can be obtained using inverse trans- 
form sampling, 

t = -1/Fo log(r) (A4) 

with r being a uniformly distribiited random number in 
the interval (0,1]. For time dependent r{t) the analytic 
solution for the probability distribution function is 

Pit) = T{t) exp ( - T{ti) dii) . (A5) 

Non-trivial time dependence of F(i) usually requires a 
numerical approach to sample t from Eqn. (A5). One 
straight-forward method would be to discretize time into 
small steps, calculate r{t) and its contribution to P{t) for 
every step, thereby numerically integrating the function 
P(t) up to a randomly chosen trigger value, at which 
point the event takes place. 

However, another method is employed here, which 
proves to be much faster and does not suffer from dis- 
cretization errors. It works for T(t) that have an upper 
bound Tm = sup(F(t)), or where such a condition can be 
enforced by introducing a cutoff. In the specific case of 
the trapped ion, all factors in Eqn. (Al) are easily limited 
by considering the energies Ej defining the trajectory and 
the excess micromotion parameters and using the peak 
neutral density. This upper bound can be adjusted after 
each collision event having affected Ej and n(x). 

The algorithm works by advancing the system by a 
time t according to Eqn. (A4) with Fq = F^. Then the 
rescaled rate 

7(t) = F(i)/F„ (A6) 

is calculated for the resulting state of the system after 

the time t and an event takes place if jit) > r, with r 
being another uniformly distributed random number in 
the interval [0, 1). If the event does not take place ("fit) < 
r) the algorithm iteratively loops back to advance the 
system by an additional t, again according to Eqn. (A4) . 
The method is exact in that it reproduces the probability 
distribution function Eqn. (A5). A proof follows below. 
The efficiency of the method is the ratio of the average of 
T{t) and F„, e = (r(t))/r„ = (7(i)). This means that 
in order to generate N events it can be expected that 
r{t) needs to be evaluated N/e times. 

Proof: We start from writing an expression for the 
probability distribution function Pg (t) obtained with the 
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suggested method. Since the final t can be the result 
of any number of iterations, -P,(t) is a sum of all these 
possibilities, Ps{t) — Ps,i + Ps,2 + Ps,3 + • • •, where P^^j 
is the probability that t results as the time of event after 
i iterations. The first few terms are given below. 

Ps,i ^ lit) PrAt) (A7) 

Ps,2=7{t) I PrMPTAt-ti){l-7{tl))dh (A8) 
Jo 

Ps,3=l{t) I r PrMPvAti-h) 
Jo Jo 

Pr„ {t-ti){l- 7(ti)) (1 - 7(f2)) dt2 dti (A9) 

Note that the products of Pr„ (as defined in Eqn. (A3)) 
in the integrals always combine to T^^^ PT^{t). There- 
fore Pr„ (t) is taken out of the sum as a common prefac- 
tor, 

Ps{t) = 7(i) Pr^ (t) (l + ^n^ ^ (1 - 7(il)) dii 

+ J / ' (1 - (1 - 7(i2)) dt2 dii + . . . ) . 

(AlO) 

Now the upper boundaries of all the partial integrals can 
be set equal to t, since they only induce ordering to the 
time series {ti,t2, ■ ■ ■ ,tn}. This rescales the terms by 
the number of possible orderings (n!). Then the partial 
integrals reduce to a single integral to the power of n, 

f r ... r""(i-7(i,))(i-7(i,))... 

Jo Jo Jo 

... (1 - 7(t„)) dtn... dt2 dh 

= ^ f /*■•• Al-7(ii)) (l-7(i2))... 

Jo Jo Jo 

...(l-7(i„))dt„... dtadii 

^^(/^'"^^'^^^^'O"- ^^^^^ 

Combining this with Eqn. (AlO) gives 

P,(i)=7(t)^r„(t)E^( / (l-7(ti))dii) 

n=0 ^-^0 ^ 

= 7(*) Pr^ (t) exp (r„,, ^ (1 - 7(*i)) dii) , (A12) 
and it follows with Eqn. (A3) and Eqn. (A6) 

P,{t) = 7(i) Tm exp (^-Tmt + TmJ^ (1 - 7(ii)) dii) 

= r(<)exp(-^ r(ti)dii) . (A13) 

This is identical to Eqn. (A5) which proves that the 
method reproduces the exact probability distribution. 



Appendix B: Structure of the simulation loop 

The following is a short description of how the main 
simulation loop has been implemented. 

1) setup system configuration: ojj , fir 
and excess micromotion, /3, neutral trap 
frequencies , Ta, Na, n{x) , ... 

2) select initial state {Ej,(pj)j^^^y z} 

3) calculate maximal collision rate T^n, 
considering nmax{Ta,Na) and Ej 

4) increment time t by Eqn. (A4) 

5) calculate ion position cuid velocity, 
using Eqn . (4) , (5) , (6) , {Ej , ipj , t) (fio„ , Vion) 

6) calculate ^{t) with Eqn. (A6), 
proceed to 7) with probability 7(t), 
else go back to 4) 

7) choose collision parameters {9, 4>) 
according to ^ 

8) update neutral atom parameters Ta, Na, 
using Eqn. (15), (16) and Etd,a 

9) apply scattering Eqn. (7) to Vion 

10) calculate new trajectory parameters, 
using Eqn. (4) , (5) , (6) , {fion,Vion, t) {Ej,<fij) 

11) loop back to 3) 

The setup of the system configuration in 1) contains 

mainly parameters which do not change during the col- 
lisions, such as trap frequencies or the atom-ion mass 
ratio. Exceptions are the neutral atom number Na and 
temperature Ta, which act as initial conditions. The ini- 
tial state for the ion energy in 2) is mostly unimportant, 
as the simulation will iterate over many collisions and 
the information of the initial state is lost after a few col- 
lisions. When looking at steady state statistics of the ion 
energy, the values after the first few collisions can sim- 
ply be ignored, thus effectively letting the system evolve 
for a short time before measuring its properties. The 
points 3) to 7) implement the collision time sampling 
algorithm described in Appendix A. The maximally pos- 
sible collision rate is calculated from the peak density 
nmax{Ta, Na) of the neutral atoms and on the maximum 
value of a{Ec) Vion{t). The latter is limited by the high- 
est possible ion velocity, which depends on the secular 
energy Ej and the excess micromotion. The calculation 
of 7(t) in 6) uses the Vion and r,on obtained in 5) to 
determine T(t) with Eqn. (9). If a collision takes place 
at the chosen time t, 7) gives the scattering angles ac- 
cording to the differential cross-section. 8) simulates the 
back-action on the neutral atoms. 9) modifies the veloc- 
ity of the ion. In 10) the new trajectory parameters Ej 
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and are determined, to represent the motion of the 
ion up to the next coUision. 

Information on any system parameter can be retrieved 
at user-defined points within the simulation loop. The 



ion energy is typically registered after 10), the collision 
and transferred energies, and Et after 8). Random 
sampling is used in 4) and 6), to determine the time of 
collision, and in 7) , for the scattering angles. 
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